%%%% Solution Algorithm for the Representative-Agent (with Exogenous Labor) Equivalent for the Aiyagari (1994,QJE)                  
%%%% Written by Orhan Torul & Oguz Oztunali
%%%% October 2017, Bogazici University, Istanbul

%%
%%
%% 0. housekeeping
clc;                                % clear screen
clear;                              % clear memory
close all;                          % close open windows
%%
%% 1. parameters and functional forms

% parameters
alpha  = 0.56;                      % capital income share from Penn World Table 9.0 (PWT)
b      = 0;                         % borrowing constraint as in Aiyagari (1994, QJE)
beta   = 89/100;                    % subjective discount factor from PWT
delta  = 5.5/100;                   % depreciation rate of physical capital from PWT
gamma  = 3/2;                       % inverse elasticity of intertemporal substitution from macroeconomics literature
varphi = 2/3;                       % Frisch elasticity of labor supply from macroeconomics literature
zbar   = 0.783472854785600;         % average effective labor productivity from the heterogeneous-agent model
hss    = 0.840226319492982;         % average hours worked at the stationary equilibrium from the heterogeneous-agent model 

% utility and (inverse) marginal utility functions
% contemporaneos utility is defined as u(c)-v(x) with the functinal forms below

u     = @(c) (c.^(1-gamma))/(1-gamma);          % utility from consumption
v     = @(h) (h.^(1+1/varphi))/(1+1/varphi);    % disutility from labor supply
up    = @(c) c.^(-gamma);                       % marginal utility of consumption
invup = @(x) x.^(-1/gamma);                     % inverse of marginal utility of consumption
vp    = @(h) h.^(1/varphi);                     % marginal disutility of labor supply
invvp = @(x) x.^(varphi);                       % inverse marginal disutility of labor supply 



%% 2. deterministic steady-state
 
% initial guess
y0 = [10,0.5,1,2,0.12]; 
% assignment of endogenous variables are as follows:
% y(1) capital (k)
% y(2) consumption (c)
% y(3) wage (w)
% y(4) real interest rate (r)

% option setting for the fsolve function
options= optimset('maxiter',999999999999','tolfun',10e-8,'tolx',10e-8);
[SS, fval,exitflag,output]=fsolve(@(y) steady_exog(y,alpha,beta,delta,gamma,varphi,zbar,hss),y0,options);

kss=SS(1);
css=SS(2);
wss=SS(3);
rss=SS(4);
lss=hss*zbar;
gdpss=kss^alpha*lss^(1-alpha);
incomess=gdpss-delta*kss;


%% 3. reporting results

disp ('steady-state for the decentralized economy:')
disp('-------------------------------------------------')
fprintf('output (y):              %.4f \n',gdpss);
fprintf('capital (k):             %.4f \n',kss);
fprintf('labor (h):               %.4f \n',hss);
fprintf('effective labor (l):     %.4f \n',lss);
fprintf('consumption (c):         %.4f \n',css);
fprintf('wage (w):                %.4f \n',wss);
fprintf('real interest rate (r):  %.4f \n',rss);
disp('-------------------------------------------------')
fprintf('utility from c:          %.4f \n',u(css));
fprintf('disutility from n:       %.4f \n',v(hss));
fprintf('contemporaneous utility: %.4f \n',u(css)-v(hss));

meanmatrix=[kss hss zbar lss gdpss css rss wss rss*kss wss*lss ];

% exporting steady-state values to LaTeX
Det_TableStationaryMeans=[meanmatrix];
columnlabelStationaryMeans = {'$K$', '$H$' ,'$Z$' ,'$L$', '$Y$', '$C$', '$r$', '$w$' ,'$rK$', '$wL$'};
matrix2latex(Det_TableStationaryMeans, 'Det_TableStationaryMeans_exo.tex', 'columnLabels', columnlabelStationaryMeans,'format', '%-6.3f', 'size', 'footnotesize');

%% 4. welfare calculations

% consumption equivalent gain calculations

% contemparenous consumption equivalent gain calculations

% assuming hours worked stays the same as in the steady-state

utilityaiyagari= -1.996968282436092;            %-1.996968282436092 refers to the average utilitarian contemporaneous utility in the heterogeneous-agent environment
utilityaiyagariss= utilityaiyagari/(1-beta);    % infinite summation of utilityaiyagari with subjective discounting     

copt=1.34;       % initial consumption guess, which is set above the required level
conv=1;          % initiating convergence
crit=0.000001;   % convergence criteria
upd = 0.0000001; % increment updating consumption guess

% indefinite forgone consumption
while conv>crit
utilityopt=(u(copt)-v(hss))/(1-beta);
conv=abs((utilityopt-utilityaiyagariss));
copt=copt-upd;
end

cfinal=copt;                                % necessary consumption level for indifference
cfinalpercent=cfinal/css;                   % necessary consumption (as a percentage of steady-state consumption) for indifference
cfinalforgonepercent= 1-cfinalpercent;      % forgone consumption (as a percentage of steady-state consumption) for indifference

% one-period forgone consumption

copt=0.148;       % initial consumption guess, which is set above the required level
conv=1;           % initiating convergence

while conv>crit
utilityopt=(u(copt)-v(hss))+beta*(u(css)-v(hss))/(1-beta);
conv=abs(utilityopt-utilityaiyagariss);
copt=copt-upd;
end

cfinalop=copt;                              % necessary consumption level for indifference
cfinaloppercent=cfinalop/css;               % necessary consumption (as a percentage of steady-state consumption) for indifference
cfinalopforgonepercent= 1-cfinaloppercent;  % forgone consumption (as a percentage of steady-state consumption) for indifference